use strict;

package Math::Polynomial::JenkinsTraub;
use Carp;
use List::Util qw(min max);

# the following statements set machine constants used
# in various parts of the program. the meaning of the
# four constants are...
# EPSILON the maximum relative representation error
#         which can be described as the smallest
#         positive floating point number such that
#         1.d0+eta is greater than 1.
# FLTMAX  the largest floating-point number.
# FLTMIN  the smallest positive floating-point number
#         if the exponent range differs in single and
#         double precision then FLTMIN and FLTMAX
#         should indicate the smaller range.
# BASE    the base of the floating-point number
#         system used.
# the values below correspond to the burroughs b6700
##	$EPSILON = .5 * FLT_BASE**(1-26);
##	$FLTMAX = 4.3e68;
##	$FLTMIN = 1.0e-45;

# ARES and MRES refer to the unit error in + and *
# respectively. they are assumed to be the same as
# eta.
use vars qw($EPSILON $ARES $MRES $FLTMIN $FLTMAX $LOW);

sub FLT_BASE() { 2 }	# shouldn't really matter?
sub COSR() { cos(94 * atan2(1,1) / 45) }
sub SINR() { sin(94 * atan2(1,1) / 45) }
sub RELAX() { .5 }
sub FXS_MUL() { 20 }

my $USE_COMPLEX = 1;

sub import {
	my $pkg = shift;
	for (@_) {
		if ($_ eq 'poly_roots_jt') {
			my $caller = caller();
			no strict 'refs';
			*{$caller.'::poly_roots_jt'} = \&poly_roots_jt;
			next;
		}
		if ($_ eq ':no_complex') {
			$USE_COMPLEX = 0;
			next;
		}
		die "Unknown package option: $_";
	}
	if ($USE_COMPLEX) {
		eval 'use Math::Complex';
	}
}

BEGIN {
	eval << 'POSIX';
use POSIX;
$EPSILON= POSIX::DBL_EPSILON;
$FLTMAX	= POSIX::DBL_MAX;
$FLTMIN = POSIX::DBL_MIN;
POSIX

	if ($@) {
		warn "POSIX not loaded: $@" if $^W;

		my $x = 1;
		while (1 + $x > 1) {
			$EPSILON = $x;
			$x *= .5;
		}

		$x = 2;
		while ($x * $x < 'Inf') {
			$FLTMAX = $x;
			$x *= $x;
		}

		$x = $FLTMAX;
		while ($x * 2 < 'Inf') {
			$FLTMAX = $x;
			$x *= 2;
		}

		$x = 1;
		while ($x) {
			$FLTMIN = $x;
			$x /= 2;
		}
		if ($FLTMIN < 1/$FLTMAX) {
			$FLTMIN = 1/$FLTMAX
		}
	}

	$MRES = $ARES = $EPSILON;
	$LOW  = $FLTMIN / $EPSILON;
}

{

my (	$a, $b, $c, $d, $e, $f, $g, $h, $a1, $a3, $a7,
	@svk, @qp, @p, @k, @qk,
	$u, $v, $sr, $si, $sz, $lz);

# variable-shift h polynomial iteration for a real zero.
# sss   - starting iterate
# nz    - number of zero found
# iflag - flag to indicate a pair of zeros near real axis.
sub real_iter {
	my ($sss) = @_;
	my ($t, $s, $omp);

	$s = $sss;
	my $j = 0;
	# main loop
	while (1) {
		# evaluate p at s
		my $pv = $qp[0] = $p[0];
		for (my $i = 1; $i < @p; $i++) {
			$qp[$i] = $pv = $pv * $s + $p[$i];
		}
		my $mp = abs($pv);

		# compute a rigorous bound on the error in evaluating p
		my $ms = abs($s);
		my $ee = ($MRES/($ARES+$MRES))*abs($qp[0]);

		for (my $i = 1; $i < @qp; $i++) {
			$ee = $ee * $ms + abs($qp[$i]);
		}

		# iteration has converged sufficiently if the
		# polynomial value is less than 20 times this bound
		if ($mp <= 20 * (($ARES+$MRES)*$ee - $MRES * $mp)) {
			$sz = [$s, 0];
			return $sss, 1, 0;
		}
		# stop iteration after 10 steps
		return $sss, 0, 0	if (++$j > 10);

		if ($j >= 2 && abs($t) <= .001 * abs($s - $t) && $mp > $omp) {
			# a cluster of zeros near the real axis has been
			# encountered return with iflag set to initiate a
			# quadratic iteration
			return $s, 0, 1;
		}

		# return if the polynomial value has increased
		# significantly
		$omp = $mp;

		# compute t, the next polynomial, and the new iterate
		my $kv = $qk[0] = $k[0];
		for (my $i = 0; $i < @k; $i++) {
			$qk[$i] = $kv = $kv * $s + $k[$i];
		}

		if (abs($kv) > abs($k[-1]) *10*$EPSILON) {
			# use the scaled form of the recurrence if the value
			# of k at s is nonzero
			$t = -$pv/$kv;
			$k[0] = $qp[0];
			for (my $i = 1; $i < @k; $i++) {
				$k[$i] = $t * $qk[$i-1] + $qp[$i];
			}
		} else {
			# use unscaled form
			@k = @qk;
			$k[0] = 0;
		}

		$kv = $k[0];
		for (my $i = 1; $i < @k; $i++) {
			$kv = $kv * $s + $k[$i];
		}

		$t = (abs($kv) > abs($k[-1])*10*$EPSILON)
			? -$pv/$kv : 0;
		$s += $t;
	}
}

# synthetic division of p by the quadratic 1,u,v
# quotient in qp, remainder in a, b
sub quad_div() {
	$#qp = $#p - 2;
	$qp[0] = $b = $p[0];
	$qp[1] = $a = $p[1] - $u * $b;
	for (my $i = 2; $i < @p; $i++) {
		($a, $b) = ($qp[$i] = $p[$i] - $u * $a - $v * $b, $a)
	}
}

# this routine calculates scalar quantities used to
# compute the next k polynomial and new estimates of
# the quadratic coefficients.
# type - integer variable set here indicating how the
# calculations are normalized to avoid overflow
sub calc_scale() {
	# synthetic division of k by the quadratic 1,u,v
	#$#qk = $#k - 2;
	$qk[0] = $d = $k[0];
	$qk[1] = $c = $k[1] - $u * $d;
	for (my $i = 2; $i < @k; $i++) {
		($c, $d) = ($qk[$i] = $k[$i] - $u * $c - $v * $d, $c)
	}

	if (	abs($c) <= abs($k[-1]) *100*$EPSILON &&
		abs($d) <= abs($k[-2]) *100*$EPSILON)
	{
		# type=3 indicates the quadratic is almost a factor of k
		return 3;
	}

	if (abs($d) >= abs($c)) {
		$e  = $a / $d;
		$f  = $c / $d;
		$g  = $u * $b;
		$h  = $v * $b;
		$a1 = $b * $f - $a;
		$a3 =($a + $g)* $e + $h * ($b/$d);
		$a7 =($f + $u)* $a + $h;
		# type=2 indicates that all formulas are divided by d
		return 2;
	}

	$e  = $a / $c;
	$f  = $d / $c;
	$g  = $u * $e;
	$h  = $v * $b;
	$a1 = $b - $a * ($d/$c);
	$a3 = $a * $e + ($h/$c+$g)*$b;
	$a7 = $a + $g * $d + $h * $f;
	# type=1 indicates that all formulas are divided by c
	return 1;
}

# computes the next k polynomials using scalars
# computed in calc_scale
sub next_k() {
	my $type = calc_scale;

	# use unscaled form of the recurrence if type is 3
	if ($type == 3) {
		$k[0] = 0;
		$k[1] = 0;
		for (my $i = 2; $i < @k; $i++) {
			$k[$i] = $qk[$i-2];
		}
		return;
	}

	if (abs($a1) <= abs(($type==1)?$b:$a) *$EPSILON * 10) {
		# if a1 is nearly zero then use a special form of the
		# recurrence
		$k[0] = 0;
		$k[1] = -$a7 * $qp[0];
		for (my $i = 2; $i < @k; $i++) {
			$k[$i] = $a3 * $qk[$i-2] - $a7 * $qp[$i-1];
		}
		return;
	}

	# use scaled form of the recurrence
	$a7 /= $a1;
	$a3 /= $a1;
	$k[0] = $qp[0];
	$k[1] = $qp[1] - $a7 * $qp[0];
	for (my $i = 2; $i < @k; $i++) {
		$k[$i] = $a3 * $qk[$i-2] - $a7 * $qp[$i-1] + $qp[$i];
	}
}

# compute new estimates of the quadratic coefficients
# using the scalars computed in calc_scale.
# use formulas appropriate to setting of type.
sub newest() {
	my ($a4, $a5, $b1, $b2, $c1, $c2, $c3, $c4, $vv, $uu, $temp);
	next_k;
	my $type = calc_scale;

	if ($type == 1) {
		$a4 = $a + $u * $b + $h * $f;
		$a5 = $c + ($u + $v * $f)*$d;
	} elsif ($type == 2) {
		$a4 = ($a+$g)*$f + $h;
		$a5 = ($f+$u)*$c + $v * $d;
	}

	# evaluate new quadratic coefficients.
	if ($type != 3) {
		$b1 =  -$k[-1] / $p[-1];
		$b2 = -($k[-2] + $b1 * $p[-2]) / $p[-1];
		$c1 = $v * $b2 * $a1;
		$c2 = $b1 * $a7;
		$c3 = $b1 * $b1 * $a3;
		$c4 = $c1 - $c2 - $c3;

		if ($temp = $a5 + $b1 * $a4 - $c4) {
			return $type,
				$u - ($u*($c3 + $c2) + $v*($b1 * $a1+$b2 * $a7))/$temp,
				$v * (1 + $c4 / $temp);
		}
	}

	# if type=3 the quadratic is zeroed
	return $type, 0, 0
}

# calculate the zeros of the quadratic a*z**2+b1*z+c.
# the quadratic formula, modified to avoid
# overflow, is used to find the larger zero if the
# zeros are real and both zeros are complex.
# the smaller real zero is found directly from the
# product of the zeros c/a.
sub quad {
	my ($a, $b, $c) = @_;
	my ($d, $e);

	!$a	and return [ $b ? -$c/$b : 0, 0], [ 0, 0 ];
	!$c	and return [ 0, 0 ], [ -$b/$a, 0 ];

	# compute discriminant avoiding overflow
	$b /= 2;
	if (abs($b) >= abs($c)) {
		$e = 1 - ($a/$b)*($c/$b);
		$d = sqrt(abs $e) * abs($b);
	} else {
		$e = $c >= 0 ? $a : -$a;
		$e = $b*($b/abs($c)) - $e;
		$d = sqrt(abs $e) * sqrt(abs $c);
	}

	if ($e < 0) {
		# complex conjugate zeros
		$e = abs($d/$a);
		$d = -$b/$a;
		return [ $d, $e ], [ $d, -$e ];
	}

	# real zeros
	$d = -$d if $b >= 0;
	$e = ($d - $b) / $a;
	$d = $e ? ($c/$e)/$a : 0;
	return [ $d, 0 ], [ $e, 0 ];
}

# variable-shift k-polynomial iteration for a
# quadratic factor converges only if the zeros are
# equimodular or nearly so.
# uu,vv - coefficients of starting quadratic
# nz - number of zero found
sub quad_iter {
	my ($uu, $vv) = @_;
	my ($relstp, $tried, $vi, $ui, $omp);

	# main loop
	($u, $v) = ($uu, $vv);
	my $j = 0;
	while (1) {
		# return if roots of the quadratic are real and not
		# close to multiple or nearly equal and  of opposite
		# sign
		($sz, $lz) = quad(1, $u, $v);
		return 0	if abs(abs($sz->[0]) - abs($lz->[0]))
					> .01* abs $lz->[0];

		# evaluate polynomial by quadratic synthetic division
		quad_div;
		my $mp = abs($a - $sz->[0] * $b) + abs($sz->[0] * $b);

		# compute a rigorous  bound on the rounding error in
		# evaluting p
		my $zm = sqrt(abs($v));
		my $t = -$sz->[0] * $b;

		my $ee = 2.*abs($qp[0]);
		for (my $i = 1; $i < @qp; $i++) {
			$ee = $ee * $zm + abs($qp[$i]);
		}
		$ee = $ee * $zm + abs($a+$t);
		$ee = 	(5.*$MRES + 4.*$ARES) * $ee -
			(5.*$MRES+2.*$ARES) * (
				abs($a + $t) + abs($b) * $zm
			) + 2.*$ARES * abs($t);

		# iteration has converged sufficiently if the
		# polynomial value is less than 20 times this bound
		return 2 if $mp <= 20*$ee;

		# stop iteration after 20 steps
		return 0 if (++$j > 20);

		if ($j >= 2 && !$tried
			&& $relstp <= .01 && $mp >= $omp)
		{
			# a cluster appears to be stalling the convergence.
			# five fixed shift steps are taken with a u,v close
			# to the cluster
			$relstp = sqrt max($relstp, $EPSILON);
			$u = $u - $u * $relstp;
			$v = $v + $v * $relstp;

			quad_div, next_k for 1 .. 5;

			$tried = 1;
			$j = 0;
		}

		$omp = $mp;
		# calculate next k polynomial and new u and v
		(undef, $ui, $vi) = newest;

		# if vi is zero the iteration is not converging
		return 0 if !$vi;

		$relstp = abs(1 - $v/$vi);
		($u, $v) = ($ui, $vi);
	}
}

# computes up to  l2  fixed shift k-polynomials,
# testing for convergence in the linear or quadratic
# case. initiates one of the variable shift
# iterations and returns with the number of zeros
# found.
# l2 - limit of fixed shift steps
# nz - number of zeros found
sub fixed_shift {
	my ($l2) = @_;
	my (	$vtry, $betas, $otv, $ss, $vi, $vv, $tss, $tvv, $ts, $s,
		$stry, $ots, $ui, $type, $svu, $ovv, $tv, $vpass, $iflag,
		$i, $svv, $oss, $betav, $nz, $spass);

	$nz = 0;
	$betav = $betas = RELAX;
	$oss = $sr;
	$ovv = $v;
	# evaluate polynomial by synthetic division
	quad_div;

	for (my $j = 0; $j < $l2; $j++) {
		($ovv, $oss) = ($vv, $ss);
		($otv, $ots) = ($tv, $ts);

		# calculate next k polynomial and estimate v
		($type, $ui, $vi) = newest;
		$vv = $vi;
		# estimate s
		$ss = $k[-1] ? -$p[-1] / $k[-1]	: 0;
		($tv, $ts) = (1, 1);

		next if !$j || $type == 3;

		# compute relative measures of convergence of s and v
		# sequences
		$tv = abs(($vv-$ovv)/$vv)	if $vv;
		$ts = abs(($ss-$oss)/$ss)	if $ss;

		# if decreasing, multiply two most recent
		# convergence measures
		$tvv = ($tv < $otv) ? $tv * $otv : 1;
		$tss = ($ts < $ots) ? $ts * $ots : 1;

		# compare with convergence criteria

		$vpass = $tvv < $betav;
		$spass = $tss < $betas;

		next unless $spass || $vpass;
		# at least one sequence has passed the convergence
		# test. store variables before iterating

		($svu, $svv, $s) = ($u, $v, $ss);
		@svk = @k;

		# choose iteration according to the fastest
		# converging sequence
		$vtry = $stry = 0;

		goto ITER_REAL	if ($spass && (!$vpass || $tss < $tvv));

	ITER_QUAD:
		$nz = quad_iter($ui, $vi);
		return $nz if $nz;

		# quadratic iteration has failed. flag that it has
		# been tried and decrease the convergence criterion.
		$vtry = 1;
		$betav *= RELAX;

		# try linear iteration if it has not been tried and
		# the s sequence is converging
		goto ITER_NEXT	if ($stry || !$spass);

		@k = @svk;
	ITER_REAL:
		($s, $nz, $iflag) = real_iter($s);
		return $nz	if $nz;

		# linear iteration has failed. flag that it has been
		# tried and decrease the convergence criterion
		$stry = 1;
		$betas *= RELAX;

		if ($iflag) {
			# if linear iteration signals an almost double real
			# zero attempt quadratic interation
			$ui = -($s+$s);
			$vi = $s * $s;
			goto ITER_QUAD;
		}

	ITER_NEXT:
		# restore variables
		($u, $v) = ($svu, $svv);
		@k = @svk;

		# try quadratic iteration if it has not been tried
		# and the v sequence is converging
		goto ITER_QUAD	if ($vpass && (!$vtry));

		# recompute qp and scalar values to continue the
		# second stage
		quad_div;
	}
	return $nz;
}

# finds the zeros of a real polynomial
# op  - double precision vector of coefficients in
#       order of decreasing powers.
# zero  - output double precision vectors of
#                real and imaginary parts of the
#                zeros.
# to change the size of polynomials which can be
# solved, reset the dimensions of the arrays in the
# common area and in the following declarations.
# the subroutine uses single precision calculations
# for scaling, bounds and error calculations. all
# calculations for the iterations are done in double
# precision.
sub rpoly {
	my (	@pt, @roots, $sc, $xxx, $max, $dx, $xx, $cc,
		$df, $yy, $bb, $aa, $t, $ff, $cnt, $bnd, @temp,
		$nm1, $x, $zerok, $factor, $i, $xm, $min, $l, $nz);
	@p = @_;
	@k = (0) x (@p - 1);
	
	# initialization of constants for shift rotation
	$yy = -($xx = sqrt(.5));
	$yy = -$xx;

	while (!$p[0]) {
		shift @p;
		@p or return [];
	}

	# remove the zeros at the origin if any
	while (!$p[-1]) {
		push @roots, [ 0, 0 ];
		--$#p;
		--$#k;
	}


NEXT_ROOT:
	my $nn = @p;
	my $n = $nn - 1;
	# start the algorithm for one zero
	if ($n <= 2) {
		return \@roots	if ($n < 1);

		# calculate the final zero or pair of zeros
		if ($n != 2) {
			push @roots, [ -$p[1] / $p[0], 0 ];
		} else {
			push @roots, quad(@p[0,1,2]);
		}
		return \@roots;
	}

	# find largest and smallest moduli of coefficients.
	($max, $min) = (0, $FLTMAX);
	for (@p) {
		$x = abs;
		$max = max($x, $max);
		$min = min($x, $min)	if $x;
	}
	# scale if there are large or very small coefficients
	# computes a scale factor to multiply the
	# coefficients of the polynomial. the scaling is done
	# to avoid overflow and to avoid undetected underflow
	# interfering with the convergence criterion.
	# the factor is a power of the base
	$sc = $LOW/$min;

	my $check;
	if ($sc > 1)	 { $check = 1 }
	elsif ($max < 10){ $check = 0 }
	else {		   $check = ($FLTMAX /$sc >= $max);
			   $sc ||= $FLTMIN;
	}

	if ($check && (1 != ($factor = $sc * sqrt FLT_BASE))) {
		$_ *= $factor	for @p;
	}

	# compute lower bound on moduli of zeros.
	@pt = map abs, @p;
	$pt[-1] = -$pt[-1];

	# compute upper estimate of bound
	$x = exp( (log(-$pt[-1]) - log($pt[0])) / $n);
	if ($pt[-2]) {
		# if newton step at the origin is better, use it.
		$xm = -$pt[-1] / $pt[-2];
		$x = min($x, $xm);
	}

	# chop the interval (0,x) until ff .le. 0
	while (1) {
		$xm = $x * .1;
		$ff = $pt[0];
		for (my $i = 1; $i < @pt; $i++) {
			$ff = $ff * $xm + $pt[$i]
		}
		last if ($ff <= 0);
		$x = $xm;
	}
	$dx = $x;

	# do newton iteration until x converges to two
	# decimal places
	until (abs($dx/$x) <= .005) {
		$ff = $pt[0];
		$df = $ff;
		for (my $i = 1; $i < $n; $i++) {
			$ff = $ff * $x + $pt[$i];
			$df = $df * $x + $ff;
		}
		$ff = $ff * $x + $pt[-1];
		$dx = $ff/$df;
		$x -= $dx;
	}
	$bnd = $x;

	# compute the derivative as the initial k polynomial
	# and do 5 steps with no shift
	$nm1 = $n - 1;
	$k[0] = $p[0];
	for (my $i = 1; $i < $n; $i++) {
		$k[$i] = $nn - (1 + $i) * $p[$i]/$n;
	}
	($aa, $bb) = ($p[-1], $p[-2]);

	$zerok = ($k[-1] == 0);
	for (1 .. 5) {
		$cc = $k[-1];
		if (!$zerok) {
			# use scaled form of recurrence if value of k at 0 is
			# nonzero
			$t = -$aa/$cc;
			for (my $i = 2; $i < $n; $i++) {
				my $j = $n - $i;
				$k[$j] = $t * $k[$j-1] + $p[$j];
			}
			$k[0] = $p[0];
			$zerok = ( abs($k[-1]) <= abs($bb) * $EPSILON * 10 );
		} else {
			# use unscaled form of recurrence
			for (my $i = 2; $i < $n; $i++) {
				my $j = $n - $i;
				$k[$j] = $k[$j-1];
			}
			$k[0] = 0;
			$zerok = !$k[-1];
		}
	}

	# save k for restarts with new shifts
	@temp = @k;

	# loop to select the quadratic  corresponding to each
	# new shift
	for (my $cnt = 1; $cnt <= 20; $cnt++) {
		# quadratic corresponds to a double shift to a
		# non-real point and its complex conjugate. the point
		# has modulus bnd and amplitude rotated by 94 degrees
		# from the previous shift
		($xx, $yy) = (	COSR * $xx - SINR * $yy,
				SINR * $xx + COSR * $yy );

		$sr = $bnd * $xx;
		$si = $bnd * $yy;

		$u = -2.0 * $sr;
		$v = $bnd;
		# second stage calculation, fixed quadratic
		$nz = fixed_shift(FXS_MUL * $cnt);
		if ($nz) {
			# the second stage jumps directly to one of the third
			# stage iterations and returns here if successful.
			# deflate the polynomial, store the zero or zeros and
			# return to the main algorithm.
			push @roots, $sz;
			push @roots, $lz if $nz > 1;

			# keep array lengths consistent
			$#qp = $#p - $nz;
			@p = @qp;
			$#k -= $nz;

			goto NEXT_ROOT;
		}
		# if the iteration is unsuccessful another quadratic
		# is chosen after restoring k
		@k = @temp;
	}
	# return with failure if no convergence with 20
	# shifts
	warn "Not all roots found" if $^W;
	return \@roots
}

}

sub poly_roots_jt(@) {
	my $r = &rpoly;
	$USE_COMPLEX
		? map+($_->[1]? Math::Complex->new(@$_) : $_->[0]), @$r
		: map+($_->[0]), grep+(!$_->[1]), @$r
}

1;

__END__
=head1 NAME

Math::Polynomial::JenkinsTraub -- Finding roots of polynomials with real
coefficients using the Jenkins-Traub algorithm.

=head1 SYNOPSIS

    use Math::Polynomial::JenkinsTraub qw(poly_roots_jt :no_complex);
    my @coef = (1, 2, 3, 1);
    print poly_roots_jt( @coef );

=head1 DESCRIPTION

This module is a translation of original Fortran code published by
M. A. Jenkins.  '$'s and '@'s are added to make it parse under Perl.

The module has one subroutine available for export: C<poly_roots_jt>.
It takes a list of numbers as the polynomial coefficients in descending
order of power, and returns a list of roots found.
Coefficents must I<not> be Math::Complex objects, not even if they
are in fact real.

There is one package option upon loading: C<:no_complex>, which tells
the module to only return real roots.  Without it, complex roots will
be returned in Math::Complex form.  Note that the algorithm still has
to find all the complex roots anyway, so this does not offer much
performance gain; rather, it is so that Math::Complex modules does not
have to be loaded if not needed.

There is an existing polynomial solving module on L<CPAN> named
Math::Polynomial::Solve.  That module offers more options in terms
of root finding; this one is generally faster.  With relatively
low degree polynomials, both are pretty accurate.  When the degree
gets high -- "high" as in "to the 10th-ish power" -- both modules
are pretty lousy.  This is due to the inherent bad conditioning
of polynomial roots and machine precision limit, so there.

=head1 REQUIRED MODULES

Standard modules List::Util and Math::Complex; POSIX is used to
obtain certain info about floating point numbers.  If you are on
a platform without POSIX, such info will be probed upon loading.

=head1 COPYRIGHT

The original Fortran code was published in
I<ACM TOMS> B<1> (1975) 178-189.
Per ACM Algorithm Policy (http://toms.acm.org/AlgPolicy.html),
ACM has the copyright on the algorithm, but "...unrestricted use
of the algorithm within a computer is permissible. General
permission to copy and distribute the algorithm without fee is
granted provided that the copies are not made or distributed for
direct commercial advantage."  Use of this perl module should
adhere to the above guideline.  There is no further restriction
by the translator.

